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We consider direct dark matter detection rates and investigate the difference between a standard 
Maxwell-Boltzmann velocity distribution and a "realistic" distribution like the ones extracted from 
numerical A''-body simulations. Sizable differences are observed when such results are compared to 
the standard Maxwell-Boltzmann distribution. For a light target both the total rate and the annual 
modulation are reduced by ~ 25%. For a heavy target the total rate is virtually unchanged, whereas 
the annual modulation is modified by up to 50%, depending on the WIMP mass and detector energy 
threshold. We also consider the effect of a possible velocity anisotropy, and the effect is found to be 
largest for a light target For the realistic velocity distribution the anisotropy may reduce the annual 
modulation, in contrast to the Maxwell-Boltzmann case. 

PACS numbers: 95.35.+d, 12.60.Jv 



INTRODUCTION 

The universe is observed to contain large amounts of dark matter [l], Q , and its contribution to 
the total energy density is estimated to be ~ 25%. This non-baryonic dark matter component, 
responsible for the growth of cosmological perturbations through gravitational instability, has still 
not been detected directly. Even though there exists firm indirect evidence from the halos of dark 
matter in galaxies and clusters of galaxies it is essential to detect matter directly. 

The possibility of direct detection, however, depends on the nature of the dark matter constituents. 
Supersymmetry naturally provides candidates for these constituents [3, 4, 5, 6]. In the most favored 
scenario of supersymmetry, the lightest supersymmetric particle (LSP) can be described as a Majo- 
rana fermion, a linear combination of the neutral components of the gauginos and higgsinos. 

Since the LSPs (or WIMPs) are expected to be very massive (towimp ~ 30 GeV) and extremely 
non-relativistic with average kinetic energy (T) ~ 50 keV (mwiMp/100 GeV), a WIMP interaction 
with a nucleus in an underground detector is not likely to produce excitation. As a result, WIMPs 
can be directly detected mainly via the recoil of a nucleus {A, Z) in elastic scattering. The event 
rate for such a process can be computed from the following ingredients: 

1. An effective Lagrangian at the elementary particle (quark) level obtained in the framework of 
the prevailing particle theory. For supersymmetry this is achieved as described in refs. @, 0]i 
for example. 

2. A well defined procedure for transforming the amplitude obtained using the previous effective 
Lagrangian from the quark to the nucleon level, i.e. a quark model for the nucleon. This step 
in SUSY models is non-trivial, since the obtained results depend crucially on the content of 
the nucleon in quarks other than u and d. 

3. Knowledge of the relevant nuclear matrix elements 0, Q, obtained with reliable many-body 
nuclear wave functions. Fortunately, in the case of the scalar coupling, which is viewed as the 
most important, the situation is a bit simpler, as only the nuclear form factor is needed. 

4. Knowledge of the WIMP density in our vicinity and its velocity distribution. Since the essential 
input here comes from the rotation curves, dark matter candidates other than the LSP are 
also characterized by similar parameters. 

In the past various velocity distributions have been considered for the dark matter gas in our 
galaxy. The most popular one is the isothermal Maxwell-Boltzmann (M-B) velocity distribution 



with {v"^) = 3t;^ ~ 3uq/2, where = {vD — {Vy) — {vD and is the velocity of the sun around 
the galaxy, i.e. vq ~ 220 km/s. Extensions of the M-B distribution have also been considered, in 
particular those that are axially symmetric with enhanced dispersion in the galactocentric direction 
[loL [ill . [H, ■ The need for such a velocity anisotropy is numerically well founded [3] ■ It has 
also been shown that the velocity anisotropy is, in principle, measurable in a direction sensitive 



experiment [15[, and it has possibly been observed to be non-zero in clusters of galaxies [16 1. 

Non-isothermal models have also been considered. Among these one should mention dark matter 
orbiting the Sun [l3| , or dark matter which is part of the Sagittarius tidal stream ^3] ■ The velocity 
distribution has also been obtained in "adiabatic" models employing the Eddington method 



velocity 

3,0, 



2l|, |22|. In such an approach, given the density of matter, one can obtain a distribution that depends 
both on the velocity and the gravitational potential. Evaluating this distribution in a given point in 
space, e.g. in our vicinity, yields the velocity distribution at that point in a self-consistent manner. 
Unfortunately this approach is applicable only if the density of matter is spherically symmetric and 
the distribution depends only on energ y. Also variants of the M-B distribution resulting from a 
coupling of dark matter to dark energy [23| have been considered. 

In the present work we will consider a Tsallis type velocity distribution 24] . This has been found to 
be a good description of the vel ocity distributions in numerical simulations with realistic dark matter 



density and anisotropy profiles (25| . We compare the direct dark matter detection rates obtained in 



this way with the results for an axially symmetric Maxwell-Boltzmann velocity distribution. 



A REALISTIC VELOCITY DISTRIBUTION 

Let us first introduce and discuss all the details of the velocity distribution function (VDF) which 
we will use. 

Very often the VDF of the dark matter particles is approximated by a Maxwell-Boltzmann (M-B) 
shape (a Gaussian). There are many good reasons for doing this, in particular the fact that many 
steps can be made analytically since the M-B is easy to integrate. It is, however, well known that 
the M-B is only an approximation, which is accurate only for an isothermal sphere. For more general 
density profiles the shape of the VDF has a different form. For example, when considering a power- 
law in density over many orders of magnitude one finds for isotropic structures that the VDF has 
the Tsallis shape [26,] . The Tsallis VDF depends on a entropic index q in such a way that the normal 



M-B is the limiting case for q ^ 1 [2J| . The shape of the VDF for more realistic cosmological dark 
matter structures has still not been calculated analytically, partly because the distribution function 
cannot be assumed to depend only on energy. First steps in the direction of actually deriving the 
VDF from a generalized Eddington method have been taken [27l |. 

It has recently been identified that the shape of the VDF actually is different for the galactic 
radial and tangential directions [25| . This means that when the particle velocity is decomposed 
into the radial and the tangential component with respect to the galactic center, the corresponding 
distributions are different. Specifically, it appears that the tangential VDF is always well fit with 
the Tsallis shape using an entropic index of g — 5/3, at least for velocities smaller than 1.6 times the 
velocity dispersion. Thus the tangential VDF is universal at any radius and the only free parameter 
is the tangential velocity dispersion. 

The radial VDF, however, differs strongly as a function of radius. It appears that the form of 
the radial VDF is fairly similar to what comes out of the Eddington method, when using a density 
profile of the NEW 28] or Sersic shape 2^, 3^] . In the inner part of the dark matter halo the shape 
of the radial VDF is fairly simple, and again the shape is reasonably well fit by the Tsallis shape 
with entropic index fairly close to unity. A clear advantage of approximating the real VDFs by 
distributions of the Tsallis shape is, as in the Maxwell-Boltzmann case, that certain integrals can 
be done analytically. 

We therefore consider a distribution function in the radial direction which is given by: 

/.(r;,a„g)=7V(9,a,) (^l-|i— , v^ = vl. (1) 



where ar is the radial velocity dispersion and N{q,(7r) is a normalization factor. We note that 
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For <7 > 1 this function is defined on the whole axis. For q < 1 the domain is limited, since one must 
demand that the function is non- negative. Thus for q = 3/4 the function becomes 
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We emphasize that this distribution thus automatically imposes an upper bound on the acceptable 
velocities. In the case of the M-B distribution this is imposed by hand, Uosc = 2.84vq, where vq is 
the solar rotational velocity. 

The corresponding Maxwell-Boltzmann distribution, i.e. the limit of ([1]) as g — s- 1, is given by: 
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Both the Tsallis shape with q ^ 3/4 and the M-B shope (corresponding to q = 1) are shown in 
Fig. 1. 

The tangential velocity distribution, defined and normalized in two dimensions, is characterized 
by g > 1 and takes the form: 
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FIG. 1: The normalized velocity distribution in the radial direction fr{v) {v = Vr) for (3 — as a 
function of vjoypl for q = 3/4 (thin line) and the one dimensional M-B distribution (thick line). 
For /3 = 0, cr = (Tr- Since this one dimensional distribution is symmetric, only the sector u > is 

shown. 



Introducing the asymmetry parameter /? defined by 



/3 = l-4, (5) 

and using the relation for the total velocity dispersion 

30-2 = al + 2a\ (6) 
we can express this distribution in terms of /? and a . Thus we obtain: 



1 1-1/5 A (9-1)1-1/3 



The corresponding two-dimensional M-B distribution with asymmetry /3 is given by: 

This can be obtained from eq. ([7]) in the limit \. As explained earlier, we assumed g = 3/4 for 
jr and g = 5/3 for /t for the more realistic VDF, which implies that the range of allowed velocities is 
V < 3a/ y/l — (2/3)/?, which is slightly different from the relation v < 2.84t;o ~ 4cr normally imposed 
on the M-B distribution. The shapes of the two 2-dimensional tangential velocity distributions are 
shown together in Fig. [2] and exhibit substantial differences in the width and the attained maximum, 
as well as at high velocities. Note that the factor 2ttv is just from the 2-dimensional phase-space 
volume. 

In terms of /3 and a the radial velocity distributions are 
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2Tra 

We note that in the case of the M-B distribution a = vo/\/2, where vq is the sun's rotational 
velocity. In the case of the Tsallis functions considered here, one has (u^) — ia^ for all q. We impose 
the condition {v^) — (3/2)wq (as for the M-B distribution) even for q 7^ 1 which implies u = wo/%/2. 

THE DIRECT DETECTION EVENT RATE 

We will now calculate both the total detection rate and the annual modulation rate. However, 
before computing the detection event rate it is necessary to transform the distributions from the 
galactic to the local coordinate system. 

From the kinematics of the WIMP-nucleus collision one has that the momentum transfer to the 
nucleus is given by 

q^2firV^, (11) 

where ^ is the cosine of the angle between the WIMP velocity and the momentum of the outgoing 
nucleus, and Hr is the reduced mass of the system. Instead of the variable ^ one can introduce the 
energy Q transferred to the nucleus, Q = q"^ / {2Amp) , where ArUp is the nuclear mass. Thus 

2{nrvy 
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FIG. 2: The 2-dimcnsional Tsallis and M-B normalized tangential velocity distributions 27ru/t(u) 
are shown as functions of v {y = Vt), for values of the asymmetry parameter P = (thick 
line) and P = 0.3 (fine line). The Tsallis curves with q = 5/3 are the ones reaching the highest 

maximum value at lower velocity. 



Furthermore, for a given energy transfer the velocity v is constrained to be 



_ j QAmp 1 
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We will find it convenient to introduce, instead of the energy transfer, the dimensionless quantity u 

. = 1(,6)^.|- , Qo = ^^4.1xlO^A-/3keV, (13) 

where b is the nuclear (harmonic oscillator) size parameter. 

It is therefore clear that for a given energy transfer the velocity is restricted from below, and 
we have already mentioned that the velocity is bounded from above by the escape velocity. We 
introduce a normalized velocity by the dimensionless variable 

y = v/ (crV2), a^/^t <y< ?/osc, (14) 

with a = {2firba) ^ . Thus, 

= ^dy. (15) 

y 



The event rate for the coherent WIMP-nucleus elastic scattering is given by [3lj : 



R = V^/coh(AMr-(^))^p\o, (16) 



or, using typical numerical values. 



i? ~ 1.60 10-3-^ ^ ^ 1 ; , /coh A . 17 

ly0.3GeVcm~3 1kg 280kms"^ 10-6 pb'' " ^ ^ 



Here, is the WIMP-nucleon scalar cross section, p(0) is the WIMP density in our vicinity, m^o 
is the WIMP mass, m is the target mass, A is the number of nucleons in the nucleus, (w^) — 3uq/2 
is the root-mean-square WIMP velocity. The dimensionless quantity /coh(^, Mr(^)) is given by 
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The quantity of interest to us is r = <coh (1 + /icohcosa), which is also dimensionless. The angle 
a describes the position of the Earth in its solar orbit. This can be viewed as a "relative" rate 
and contains all the information regarding the WIMP velocity distribution and the structure of the 
nucleus. It also depends on the reduced mass of the system. 

The event rate is prop ortional to the WIMP flux, i.e. proportional to the WIMP velocity. It is 
not difficult to show [SlJ that 

^ = F'{u) .^y^y^dy [ de f dc^friy, d) Mv , e, <P, a)e(yL - A), (19) 

du Ja^ V 3 2/^ Jq Jo 

with F(u) the nuclear form factor. In the integrand we have explicitly displayed all the factors of 
y in order to keep track of their origin. The first one comes from the WIMP flux, the second from 
the transformation and the last is the usual phase-space factor. 8(a;) is the Heavyside function 
introduced to guarantee that the WIMP velocity is bounded by the escape velocity Thus, 

/fc = + ?/g + (5o + yo^o cos a + 2ySo cos 9 sin a 

Even though the values of yesc are different, the effect of this constraint is small, with the possible 
exception of the event rates for very high threshold energy. In that special case the difference in the 
upper cutoff, ^csc, between the two models may yield significant differences. 
This way the differential event rate in eq. (|19p can be cast in the form: 



(20) 



We have seen that the parameter a depends on the nucleus, the velocity distribution [vq or a) and 
the WIMP mass. 

By performing a Fourier analysis of the function "^{xja) which is a periodic function of a and 
keeping the dominant terms we find: 



dr PI 

— = \ -a^F'^{u)-^o[a^/u) \l + H [a^Ju) cos a] 
du \ 3 



(21) 



Sometimes we will consider each term separately in the above expression by writing: 



dr dt dh dt /2 n^,, , ^, dh /2 n^, , , 

^ = ^ + ^cosa, — = ^^-a'F^{u)■^o{aV^) , ^ = \j -a^F^{u)-^o{aVu)H{a^) (22) 

It is thus clear, that and t are related to the total rates, whereas H and h are related to the 
annual modulations. 

Before proceeding further by considering a special target, it is instructive to concentrate on '^o{x) 
and the modulation H{x) as functions of the asymmetry parameter (3. For this purpose we show 
the function ^'o(a;) in Fig. [3] and the function H{x) in Fig. [4l In the case of ^'o(aY^), one clearly 
sees that the effect of the asymmetry parameter /3 is small for all energy transfers. In the case of the 
M-B distribution one has that larger (3 imply smaller differential rate, while in the case of the more 
realistic distribution a larger asymmetry implies a larger differential rate. In both cases the function 
H(ay/u) is an increasing function of accompanied by a change in sign. Thus in obtaining the 
modulation of the total amplitude the contribution of the low Q section tends to cancel the high Q 
part. Which part will dominate depends on the reduced mass and the nuclear form factor. In the 
case of the present realistic distribution, the velocity asymmetry has a small effect on the modulation 
of the differential rate, and this effect occurs at high energy transfers. 



(a) (b) 
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FIG. 3: The time independent part of the differential rate, '^o{ay/u), as a function of the energy 
transfer ay/u. The normaUy used M-B distribution is on the left and the present realistic velocity 
distribution is on the right. There is essentially no dependence on the velocity asymmetry 

parameter. 
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FIG. 4: The same as in Fig. [3] for the function H{ay/u), which describes the coefficient of the time 
dependent part of the differential rate (differential modulated amplitude). Note the change in sign 

as a^/u, i.e. the energy transfer, increases. This leads to cancellations in the total modulation 
amplitude. M-B is on the left, Tsallis is on the right. In some cases it may even become negative 
(maximum event rate in December). The dependence on the asymmetry is small, visible only at 

high energy transfers. 



APPLICATIONS 



Let us now focus on the aspects affected by the WIMP velocity distribution. The total (time 
averaged) rate is given by: 



''coli 



^du, (23) 
du 



where ■Umin is determined by the detector threshold and Mmax = {VcscY /o? by the maximum WIMP 
velocity. By including both "^^[a^Ju) and H{ay/u) we can cast the rate in the form: 

'^coh = icoh (1 + hcoh cos a) 

1 /""ma" dhcoh 



du, (24) 
where the dimensionless quantity r characterizes the total rate, and h is the annual modulation. 
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FIG. 5: (a) The form factor F'^iu) for ^^F empfoyed in our calculation as a function of u = Q/Qo, 
where Q is the energy transfer to the nucleus and Qo = 809 keV. (b) The same quantity as a 

function of the energy transfer Q. 



The direct WIMP detection rate depends on the nucleus via its form factor and its mass. It also 
depends on the WIMP mass through the reduced mass fir, entering through the parameter a. For 
our numerical study we will focus on two very different targets to show general trends, namely a 
medium heavy target, ^^^I, and and a light one, ^^F, which are two of the most popular targets 
employed. 



The case of a light target 



The actual results have been obtained for ^^F, but those of other light targets are similar. The 
nuclear form factor we use was obtained in the shell model description of the target and is shown in 
fig. m Integrating over the energy transfer, assuming either no detector cut off (wmin = 0) or a cut 
off of Qthr = 5 keV, we obtain the total rate in the case of the target "'^^F. The results are shown in 
Fig. [6] for t and in Fig. [7] for h. 



The case of an intermediate mass target 

Even though the actual results have been obtained for ^^^I, the situation is similar for other 
intermediate targets. The nuclear form factor employed was obtained in the shell model description 
of the target and is shown in Fig. [8] In the case of the target ^^^I the obtained results are shown in 
Figs. [9l and fTOl for Qmin = at the top and Qniin = 10 keV at the bottom. 



DISCUSSION 



In our formalism the quantities which are affected by the velocity distribution, are the parameters 
t, describing the total time-averaged rate, and h describing the amplitude of the time dependent 
rate (modulation effect). These parameters also depend on the target nucleus as well as the WIMP 
mass. In this work they have been obtained in the case of the coherent scattering, but they are 
not expected to be very different in the case of the spin mode. Our results can be summarized as 
follows: 

1. The time average rate. 

The parameter tcoh in the case of symmetric velocity distribution (/3 = 0) is not very different 
from that obtained in various other models, like those obtained in the Eddington approach 
HHiHil. The obtained results depend on the choice of the target. Thus: 



The parameter t for a light target 

It is at first a decreasing function of the WIMP mass and eventually becomes a constant 
for heavy WIMPs, which for the realistic VDF is approximately 20 percent lower than 
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FIG. 6: The dimensionless quantity icoh in the case of the target "'^^F for Qthr = at the top and 
Qmin = 5 keV at the bottom. The abscissa is for different WIMP masses, m^. icoh represents the 

time-independent part of the detection rate. On the left the results obtained for a M-B 
distribution, while on the right those for the more realistic distribution. The thick and thin curves 

correspond to /3 = and 0.3 respectively. 



for the M-B case. This limiting value is not much affected by threshold effects. The value 
of the parameter t at low WIMP masses is, as expected, sensitive to the attained energy 
threshold. The dependence on the asymmetry parameter is small (few percent) in both 
distributions. 

• The parameter t for a heavy target. 

The general behavior is similar to the above except that now the limiting value is an 
order of magnitude lower. This may partly offset the advantages offered by the favorable 
nuclear mass dependence of the event rate, which is proportional to ~ ^^A « (see 
Eqs. p6|) - p8)) ). not included in the plots. A threshold of 10 keV further reduces the rate 
by about a factor of two. At low WIMP mass the advantage of employing an intermediate 
target is expected, since t is only about a factor of three down from that of a light target. 
The effect of the realistic VDF is small (from few to 10 percent). The effect of asymmetry 
is close to unobservable. 

2. The modulation effect. 

We can now focus on two issues: 

• The modulation of the differential rate. 

All the previous velocity distributions imply a modulation amplitude which i) increases 
with energy transfer, and ii) changes sign as we go from low to high energy transfers. 
Also, iii) the slope and the point of the change of sign depend on the WIMP mass and 
the nature of the target. 

In the case of the present realistic velocity all the above hold true. The effect of the 
velocity asymmetry is quite small. 
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FIG. 7: The same as in Fig. [5] for the quantity hcoh, which represents the time-dependent part of 

the detection rate (the annual modulation). 




(a) (b) 

FIG. 8: (a) The form factor for ^^''l employed in our calculation as a function of m = Q/Qo, 

where Q is the energy transfer to the nucleus and Qo = 64 keV. (b) The same quantity as a 

function of the energy transfer Q. 



• The total modulation amplitude: 

i) In the case of a light target both distributions predict a similar modulation amplitude. 
Its sign is positive, i.e. the maximum is around June 3rd. The effect of asymmetry is small, 
but, when comparing the two models, it tends to go in opposite directions. Unfortunately 
the predicted difference between the maximum and the minimum at zero threshold cannot 
exceed 4%. It may double for an energy threshold of about 5 keV. 

ii) In the case of a heavy-intermediate target the modulation becomes negative for both 
models as the WIMP mass increases. Thus the event rate for both models may attain a 
maximum in December. The difference between the maximum and the minimum is now 
quite a bit higher, i.e. it can be as high as 9% depending on the WIMP mass. For a light 
WIMP mass the modulation is quite small at zero threshold. In the presence of energy 
thresholds it tends to increase. This is understandable, since the modulation here is the 




FIG. 10: The same as in Fig. [5] in the case of ftcoh- The dependence on the asymmetry parameter 

fi is barely visible. 



ratio of two quantities (see Eg. [24]) ) and the denominator, which determines the average 
event rate, falls much faster then the numerator in the presence of a threshold. The 
modulation can be quite sizable, especially for light WIMPs. With a 10 keV threshold we 
predict differences between maximum and minimum of about 5-15% in the case of both 
an M-B distribution and the present distribution. The event rate may attain a maximum 
in June for light WIMP's. The effect of the velocity asymmetry is virtually negligible. 



CONCLUSIONS 



The event rate for direct WIMP detection depends on the nucleon cross section (which in turn 
depends on the assumed particle model), the WIMP mass, the structure of the target, and the 
WIMP velocity distribution. In the present study we focused on the last aspect by considering a 
realistic velocity distribution similar to the ones extracted from numerical TV-body simulations, here 
parametrized by the Tsallis shape. Generally speaking the calculations using the M-B distribution 
and the present distribution yield similar results. This is particularly true of the time independent 
part of the event rate. Some differences show up, however, in the prediction of the time-dependent 
event rate, described by the modulation amplitude h (the difference between the maximum and 
the minimum event rates being 2|ft,|). For intermediate targets at zero threshold the modulation 
amplitude obtained in the present calculation is small for light WIMPs, \h\ < 2%, but it can reach 
values of \h\ « 5% for heavy WIMPs. For light WIMPs in the presence of an energy threshold 
of a few keV the obtained values of \h\ may become larger, but this occurs at the expense of the 
number of counts. The predicted effect on the event rates of the asymmetry parameter of the velocity 
distribution is relatively small. 

It will be interesting to explore the possible consequences of the realistic velocity distribution 
considered in this work on the event rates of experiments which measure the direction of the recoiling 



nucleus j32l. l33l. |34|. 
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